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A method for measuring the pair interaction potential between colloidal particles 
by extrapolation measurement of collective structure to infinite dilution is presented 
and explored using simulation and experiment. The method is particularly well 
suited to systems in which the colloid is fluorescent and refractive index matched 
with the solvent. The method involves characterizing the potential of mean force 



between colloidal particles in suspension by measurement of the radial distribution 

a 

function using 3D direct visualization. The potentials of mean force are extrapolated 

O 

to infinite dilution to yield an estimate of the pair interaction potential, U(r). We 

use Monte Carlo (MC) simulation to test and establish our methodology as well as to 
explore the effects of polydispersity on the accuracy. We use poly-12-hydroxystearic 
^ acid-stabilized poly(methyl methacrylate) (PHSA-PMM A) particles dispersed in the 

solvent dioctyl phthalate (DOP) to test the method and assess its accuracy for three 
different repulsive systems for which the range has been manipulated by addition of 
> electrolyte. 
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I. INTRODUCTION 

Colloidal systems may undergo crystallization because of the effects of packing and ex- 
cluded volume^, repulsive or attractive charge^, and/or weak attractions caused by addi- 
tives such as non-absorbing polymer 4 . Moreover, the phase behavior can often be dramat- 
ically changed by making only small modifications to the particle-particle interactions^^. 
Knowledge of the interaction between colloidal particles and the ability to tune these in- 
teractions is important for designing and assembling target materials. Target materials, 
such as self-assembled arrays of colloids, have potential applications as periodic dielectrics^, 
photonic band gap materials^', and chemical and biological sensors^. The pair interaction 
potential, U(r), which characterizes the potential energy change that results as two isolated 
particles are brought from an infinite to a finite separation, is a common parameterization of 
the interaction between colloids. The ability to parameterize particle interactions into pair 
potentials allows simulation and theory to be readily incorporated into the experimental 
design process. Simulation and theory have been shown to be valuable tools for predicting 
and explaining structures and trends in colloidal systems. For example, simulation and the- 
ory have been used to calculate the crystal-nucleation rate of hard-sphere colloids-^, phase 
behavior of attractive^ and repulsive^ colloids, and the stability of binary ionic colloidal 
crystals^. Thus, characterization of U(r) for colloidal particles is an important step in the 
process of predictably assembling target phases. 

Common techniques for measuring the interaction potential of colloidal particles include 
colloid probe atomic force microscopy (AFM) 11 , surface force apparatus^, total internal 
reflection microscopy (TIRM)P1, and optical tweezers^. Each of these techniques allows for 
the direct measurement of the pair interaction potential between isolated particles and sur- 
faces. While these techniques provide valuable information about colloidal interactions their 
applications to the problems of colloidal assembly are limited because, with the exception of 
optical tweezers* 15 * 16 !, they do not characterize particle-particle interactions directly, which is 
often more relevant to the study of self-assembly of bulk colloids. Additionally, TIRM and 
optical tweezers are not generally applicable to refractive index matched colloidal systems 
since these technique require refractive index contrast for optimal performance. Often, col- 
loidal systems that are useful for self-assembly are approximately refractive index matched, 
which minimizes strong attractive interactions due to van der Waals forces that tend to 



trigger gelation and irreversible aggregation^. 

For colloids that are refractive-index matched, fluorescent and approximately 1/im in 
size, fluorescence or confocal microscopy methods can be used to characterize structure and 
order in such suspensions^. The ability to directly quantify the structure enables the use of 
statistical physics methods to characterize the interaction potential. This technique involves 
calculating the radial distribution function, g(r), and extracting the potential of mean force, 
W(r), using the following relationship^, 

W(r)/k b T = -\n[g(r)} (1) 

where k^ is the Boltzmann constant and T is the temperature. W(r) is a volume fraction 
(i.e. density) dependent measure of interaction; to determine the pair potential U(r), and 
not simply a potential of mean force, an extrapolation to the limit of infinite dilution is 
required, as we will discuss in this work. These methods are complementary to the pre- 
viously discussed ones because they can often be performed on the exact system that will 
be used in self-assembly (i.e. a bulk solution of colloids). This basic treatment appears in 
various implementations in the literature. Several groups have examined bulk solutions of 
colloids, capturing particle positions using video microscopy^ and confocal laser scanning 
microscopy^^. In general, these groups calculate g(r) from the microscopy data of a sus- 
pension at a particular "dilute" concentration and apply equation [T]to arrive at an "effective" 
pair potential at a finite volume fractiorpSEu That is, the method of these papers assumes 
that at dilute concentrations the potential of mean force approximates the pair potential. 
Hsu et alP 21 approached the problem by using bright-field microscopy. Using a quasi-2D 
methodology, they captured many statistically independent images of particles interacting 
with one another, and then computed the 2D radial distribution function to extract the pair 
potential. To corroborate their results, they employed Monte Carlo (MC) simulations to 
calculate g(r). Wu and BevarP^took advantage of TIRM and video microscopy to capture 
the interaction forces between particles. In this case, colloidal pair interactions were assessed 
by separating out the additional contribution of the surface that is present in TIRM studies. 
Finally, Lu et al. 24 , to support studies of gelation, parameterized short-range attractive po- 
tentials of micron-sized colloids by comparing measurements of the second virial coefficient 
and the cluster mass distributions to the results of MC and molecular dynamics simulations. 
Overall, this class of techniques has been shown to be successful in many applications, pro- 



viding results that are both consistent with theoretical behavior and simulation results. A 
drawback of this basic methodology is that there is no clear definition of what is "dilute" 
and, as we discuss in detail in this paper, if we do not sample in the correct regime this 
method can produce results that appear qualitatively correct but are quantitatively wrong. 
Moreover, this paper shows that the quantitative effect of the extrapolation methodology 
we propose should not be underestimated and that measurements of collective structure 
converted to effective potentials at a small but still finite volume fraction, i.e. assuming 
W(r) = U(r) at low volume fraction, are prone to significant systematic error. 

In this work, we outline a general procedure for the determination of the pair interaction 
potential, U(r). In this procedure we use linear regression to extrapolate the potentials of 
mean force, W(r), at finite concentrations to infinite dilution. We use both simulation and 
experiment to assess the validity of this method and provide guidelines for its use. In section 
II A we detail our method for determining the pair potential, providing a rational basis from 



statistical mechanics. In section II B| we introduce our simulation model and method, and 



in section II C we introduce our experimental method. In section HI A we use simulation 



to determine the dilute limit, i.e. the regime where the two body forces are dominant and 



W(r) scales linearly with volume fraction. In section III B we use simulation to test the 



accuracy of the extrapolation method. In section |IHC| we explore the role of polydispersity 
as it affects the potential derivation. In section IV we apply the guidelines from section 
II to an experimental system of fluorescently labeled poly-methyl methacrylate (PMMA) 
particles stabilized by poly-12-hydroxystearic acid (PHSA) and compare to the theoretical 
screened Coulombic potential. In Section M we provide concluding remarks. 



II. METHODOLOGIES 

A. Method for determining U(r) from the potentials of mean force 

Following Chandler^, we can develop a rational basis for our method from statistical 
mechanics. Starting with: 

9(r) = e-^M (2) 

we have a relationship between the radial distribution function, g(r), and the potential 
of mean force, W(r). W(r) represents the reversible work for the process of moving two 



particles from infinite separation to a finite separation of r. W(r) can be separated into two 
parts: 

W{r) = U(r) + AW(r) (3) 

where U(r), the pair potential, is the reversible work to move two isolated particles to a 
separation of r at infinite dilution, and AW(r) is the contribution to W(r) due to the 
density of the system (i.e. interactions with surrounding particles in the system, not many 
body forces. This can additionally be thought of as moving through a potential landscape 
resulting from the presence of other particles). Combining equations |2J and [3] and taking the 
natural log, we arrive at: 

- ln i^ r) = U{r) + AW(r) (4) 

In the limit of infinite dilution (i.e. when the volume fraction, <p, goes to zero) we have 
lim^o AW(r) = 0, and thus: 

lim ~ ln f^ = lim W(r) = U{r) (5) 

</>-»0 p 0->O 

Equation [5] is the basis for our method for determining U(r). To avoid confusion, we 
will now refer to the potential of mean force as W(r, <ft) and the radial distribution function 
as g(r,(j)), to highlight the volume fraction dependence. The general procedure to calculate 
U{r) is as follows: 

1. Calculate g(r, 4>) at a series of finite, dilute values of (p 

2. Calculate W(r, <p) from g(r, 0), using equation [IJ for each value of <p 

3. For each value of r, perform a linear regression of W(r, 0) vs. <fi 

4. For each value of r, evaluate the linear regression at <fi = to construct an estimate 
of U(r) 

We should note that unlike other methods in the literature, this method does not rely on a 
single g(r, <fi) but rather on the behavior of a collection of g(r, <fi) data over a range of 0. In 
section II we will explore the benefits of using a collection of g{r,<f)) data and will provide 
guidelines for the use of this method. It is important to note that this methodology only 



applies to systems where many body forces are negligible. It has been shown that when 
the electrostatic screening length is less than the interparticle separation, the assumption of 
pairwise additivity is well suited to predict the forces between colloids^ the experimental 
systems we study in t his paper fit this criteria. It is additionally important to note that the 
charge of the colloids may be a function of 0, thus making the application of U(r) to higher 
density systems problematic 2 -^ however, this is an issue with all methods that determine the 
interaction at low <p and will not be specifically addressed in this article. 

B. Simulation Method and Model 

To study the various aspects of colloidal pair interactions, we performed NVT Monte 
Carlo (MC) simulations that employ the Metropolis sampling algorithm 28 . Simulation is 
powerful in this application since we explicitly know the "true" interaction potential (i.e. 
it is programmed into the simulation code), have exact control over polydispersity and 
volume fraction, and have no artifacts associated with identification of the particle centroid. 
Thus we can assess the accuracy of the potential derivation method under well controlled, 
ideal conditions. We conducted simulations of spherical particles that interact via screened 
Coulombic interactions, modeled using the Yukawa potentiaP, given as: 

U(r) exp[-Ka(r/a - 1)] 

k b T e r/a l } 

where k is the inverse Debye length, a is the particle diameter, and e is the energy at contact 
(dimensionless, scaled by kbT). In our simulations, particles were treated as hardcore and 
not allowed to overlap (i.e. U(r)/kbT = oo when i/a < 1); the potential was truncated at a 
distance r when U(r)/k(,T < e/60^. In the limit of very large k, the Yukawa potential be- 
comes very short-ranged and essentially reduces to the hard sphere potential, only capturing 
excluded volume^. We used the hard sphere potential, given by equation |7| to investigate 
this limit. 
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To model polydispersity in particle size, particle diameters were set based on a prescribed 
Gaussian distribution where the average particle diameter is given as a. Particle interactions 



were modeled by Equations [6] and [7j however the radial separation, r, was adjusted based 
upon the deviation from the mean diameter. Specifically for a pair of particles i and j, with 
diameters <Tj and <Tj, respectively, and a center-to-center separation of r, the scaled radius, 
r sca i ed , is given by: 

r scaled = r - 0.5(ct; -a)- 0.5(0, - a) (8) 

Note that eqn [8] reduces to r sca ied = t when <7; = <7j = a. Thus, the Yukawa potential for 
a polydisperse set of particles is given as: 



Uij(r) exp[-/scr([r + a- 0.5ctj- 0.5<7j-)]/cr- 1] 

where again, a is the mean particle diameter. This fixes the interaction range of the potential 
with respect to the surface of the particle. This assumption is reasonable since the Debye 
length (1/k) is the dominant parameter in determining interaction range and is independent 
of particle diameter. Additionally, for simplicity, we make the assumption that the energy at 
contact is not a function of particle diameter, instead focusing primarily on the role excluded 
volume interactions play in polydispersity. Such an assumption appears reasonable; with a 
fixed e value, a 5% increase or decrease in particle diameter only results in a ~2.5% change 
in the charge number, Z, and a 10% change only results in ~5% change in Z, as calculated 



using equation 10 



In all cases, we used system sizes of 1000 spherical particles in a periodic box and fixed 
the dimensionless temperature at a value of 1. For each system we started from a random, 
disordered configuration and allowed the system to run for approximately one million MC 
timesteps, collecting data every 1000 MC timesteps. System size effects should be minimal 
considering the low volume fractions studied and relatively short-ranged interaction; the 
box length was substantially larger than the range of interaction for all conditions studied 
(e.g. for (p = 0.1 the box length, lt, ox > 17a and for = 0.001, lb ox > 80a). Dimensionless 
temperature is defined as T* = k T /e. The interaction potential programmed into the 
simulation code will be referred to as the "known" interaction potential. 



C. Experimental Method 

1. Synthesis and Characterization of PHS A- stabilized PMMA Particles 

Non-aqueous solvents are often selected for colloidal assembly because they offer the pos- 
sibility of matching the solvent and particle refractive index, thereby minimizing attractive 
van der Waals interactions that can interfere with assembly. A common choice for assembly 
includes mixtures of cyclohexyl bromide and decalirpH An alternative choice, particularly 
useful for confocal microscopy, rheology and field-assisted assembly is dioctyl phthalat d 30 1 31 1 . 
Because of its high viscosity, this solvent is particularly compatible with the scan rate of 
confocal microscopy, and we use it to test our method for that reason. 

Fluorescently-labeled poly(methyl methacrylate) particles stabilized by poly-12-hydroxysteric 
acid (PHSA) were synthesized using an adaptation of the methods of Antl et al.^, Campbell 
and Bartlett 33 , and Pathmamanoharan et al; 34 , as discussed by Solomon and SolomorP^. 
Particles were labeled using Nile red dye. Previous work has shown that PHSA-PMMA 
colloids in cyclohexyl bromide (CHB) and decalin are charged^. Shereda et al. also reported 
particle charging for these particles in DOP^H. Such colloidal systems have been modeled 
by the Yukawa pair potential (equation [6]) where we explicitly define: 



Z 2 X 



(10) 



a p k B T(l + K(j p /2) 2 

where Z is the charge number, A# is the Bjerrum length (11 nm for the system studied 
here23), anc [ K j s ^g j nverse Debye lengttP. 

Parameters in equations [6] and 10, estimated by experiment, are reported in Table II C 1 



/t _1 was estimated from conductivity measurements. Briefly, solutions of dioctyl phthalate 
(DOP; Sigma-Aldrich, used as received) containing 10 /jM. and 2mM of tetrabutylammo- 
nium chloride (TBAC) salt were prepared and their conductivity measured with devices 
(Model 1154 Precision Conductivity Meter, Emcee Electronics, Venice, Florida or Model 
EW-01481-61, Cole-Parmer, USA) whose performance was verified by measurements with 
known standards. As per the method of reference^, using Waldens rule^"^ (viscosity of 
DOP = 71 mPa s at 25CP^ and reference ion mobilities in water at 25C (as taken from 
reference^ TBA + ion mobility is 19.4 cm 2 S mol -1 and Cl~ is 76.3 cm 2 S mol -1 ), the ion 



concentration was determined and calculated using equation 11 



k = ^MX B pi (11) 



where pi is the density of the cation or anion 29 . The Debye length of salt-free DOP was 
not determined in this way because of the uncertain identification of the mobile ions in 
this solvent. Thus, in this case, a one-parameter fit to the extrapolated pair potential (as 
discussed in the results) was made and a value of 440 nm was obtained. This value is 
included in the table for completeness. This Debye length is qualitatively consistent with 
literature reports^ of the very low conductivity (~ 2.3 x 1CT 9 S/m) of dioctyl phthalate. 

The charge number of the PHSA-stabilized PMMA particles was determined from mea- 
surements of their electrophoretic mobility (Zetasizer Nano Series, Malvern, United King- 
dom). Particles were prepared as 1 vol percent solutions in pure DOP solvent (i.e. [TBAC] 
= 0), [TBAC] = 10/iM, and [TBAC] = 2mM. Solutions were placed in a dip cell designed for 
non-aqueous solvents. A voltage of 50mV was applied to each sample and three independent 
samples were studied. Using the calculated Debye length and measured mobility, the zeta 
potential was determined by the method of O'Brien and White^M The charge number on 
the particles was determined from: 

Z=g = 47re e r (q/2)C(1 + W2) (12) 

e e 

where q is the particle charge, e is the charge on an electron, eo is the permittivity in a 
vacuum, e r is the relative dielectric constant of DOP ( 5.1Cp^), and C, is the zeta potential. 
The distribution of particle diameters, determined from scanning electron microscopy 
(SEM) analysis of 250 particles, is plotted in Figure HI From the best fit of the normal dis- 
tribution, we determined the mean diameter to be 951nm ± 38nm. Because PMMA colloids 
may swell in organic solvents^, we compared the SEM diameter to direct measurements of 
the colloid size in solution. To perform the comparison, we prepared a sample at a very high 
TBAC salt concentration of 50mM to induce aggregation of the particles in the solvent. A 
CLSM image volume was acquired and the separation between particle pairs was computed 



by image processing per the method described subsequently in section |HC3| We found the 
particles swell approximately 5% to 1001 ± 30 nm. This small change in diameter does not 
have a major effect on our measurements or U(r). For consistency, we will report our find- 
ings based on the particle diameter measured by the CLSM method, since our experiments 



were done in DOP solvent. 

50 
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FIG. 1. Distribution of PHSA-PMMA particles diameters with fitted Gaussian curve shown as 
solid line; particles have a diameter of 951 ± 38nm. Inset is an SEM image of PHSA-stabilized 
PMMA particles used to generate the distribution; scale bar represents 2/iin. 
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No TBAC 


440 


-28 


-100 


27 


2.2 


10 fM TBAC 


250 


-30 


-160 


31 


4.0 


2 mM TBAC 


49 


-32 


-630 


35 


20 



TABLE I. Electrokinetic measurements for the three experimental systems studied. 



2. Sample preparation and image volume collection 

In this work, we explored three different systems where the solvent and particles are 
approximately refractive index (RI) matched. The PMMA-PHSA particles (RI = 1.489) 
were dispersed in DOP (RI = 1.485) containing no salt, 10/xM TBAC salt or 2mM TBAC 
salt. Samples were prepared at nominal volume fractions of </> = [0.005, 0.01, 0.015, 0.02, 0.03, 
0.04, 0.05]; however, precise estimation of the volume fraction was made from the results 



of quantitative image processing discussed in section |HC3| The samples were initially 

10 



mixed and then allowed to equilibrate for 24 hours. Samples were then gently remixed 
and subsequently transferred to glass specimen vials (outer diameter, O.D. = 12mm) that 
were adhered to a microscope cover slip using ultraviolet bonding glue (Dymax Corporation, 
United States). The cover slip was attached to a 35mm O.D. glass ring made from Pyrex 
standard wall tubing. To assess the stability of this colloidal system, samples prepared at <p 
= 0.005 in the DOP solvent were monitored for 24 hours for signs of phase instability. No 
aggregation or phase instability was observed. 

Sample imaging was performed on a Leica TCS SP2 confocal microscope (Leica Microsys- 
tems, Wetzlar, Germany). A lOOx oil immersion objective with numerical aperture 1.4 was 
used. The particles were dyed with Nile red (Sigma-Aldrich, United States) and were ex- 
cited with a green neon (GreNe) laser beam (Ao = 543nm). Emission from 550nm to 650nm 
was collected. To avoid possible effects of sample boundaries on particle configurations, all 
points in the image volumes were located at least 20/im from any boundary of the specimen 
vial. Stacks of 247 images with a resolution of 512 x 512 pixels were obtained and processed 
to extract particle centroids, as described in the next section. Images were acquired with a 
spatial resolution of 69.2 x 69.2 nm/pixel in the objective plane and an axial separation of 
81.4 nm. Thus, the size of the image volume was 35 x 35 x 20 fim 3 . 

3. Image Processing 

To identify particle locations, we used image processing algorithms based on the work 
of Crocker and GrierSl expanded to 3D systems^ as discussed in Dibble et alPH First, a 
Gaussian filter was applied to the 3D image volume. Second, particle centers were identified 
using a local brightness maximum criterion. That is, a voxel was identified as a candidate 
centroid if it was the brightest within a local cubic region of half-width w. For our systems 
and imaging conditions, w = 7, which corresponded to approximately 485 nm. Finally, using 
the moments of the local intensity distribution, particle positions were refined to subpixel 
accuracy^. For our setup, the accuracy was calculated to be ±35 nm in the x — y plane 
and ±45 nm in the z direction in reference^. g(r, 4>) was calculated using the centroid loca- 
tions. We assessed the accuracy of the image processing algorithm by examining composite 
images for which centroid locations were overlaid on the fluorescence images; the algorithm 
was found to identify nearly all of the particle centroids. We also validated the centroid 

11 



calculation for particles at close contact, finding no evidence of erroneous behavior, as has 
been observed for other systems^®'. For our fluorescently labeled system, CLSM provided 
a strong contrast between particle and solvent, and does not produce an Airy pattern; thus 
we meet both conditions shown to be essential for accurate centroid determination outlined 
in reference^. 



III. SIMULATION RESULTS 

A. Determining the Dilute Regime 

The general method for calculating the pair potential involves linearly regressing W(r, 0) 
at finite values of 0, and evaluating the regression at = 0. In order for this to be valid, 
we must be in a range where W(r, 0) scales linearly with (i.e. the regime where only 
2-body forces are dominant and there are no strong volume fraction effects). Unfortunately, 
there is no general rule for what value of constitutes this dilute regime; the crossover 
between "dilute" and "concentrated" will be a function of the strength and relative range 
of the interaction. To determine the crossover, 4> max , we performed a series of simulations 
using the hard-sphere and Yukawa potentials over a range of 0. Specifically for the Yukawa 
potential, we performed simulations with e = 3 and ko = (5.0, 3.0, 1.5, 1.0, 0.5, 0.25), e = 
50 and na = (5.0, 3.0, 1.5, 1.0), and e= 100 with na = (5.0 ,3.0, 1.5, 1.0) at a variety of 
values of 0. This allows us to explore the dependence of ma z on both interaction range and 
interaction strength. 

To determine max , the maximum volume fraction for which W(r, 0) scales linearly with 
0, we performed the following procedure. For each value of r/cr, we plotted — ln[g(r, 0)] (i.e. 
W(r, 0)) as a function of 0; we included the pair potential as the value that occurs at = 
0. We then determined, by eye, the approximate crossover between linear and non-linear 
behavior. We then performed a linear regression on the subset of the data that appeared 
linear and fine-tuned our results. We looked at both the residual of the fitting of the data 
and compared between the known (i.e. the potential programmed into the simulation) and 
regressed value at = to determine the appropriate value of maa; . Figure [2]^ plots an 
example of the crossover between "dilute" and "concentrated" regimes showing a region with 
distinct linear behavior; the extrapolation is also plotted showing near perfectly agreement 
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FIG. 2. For a Yukawa system with e = 3, na = 1.5, (A) determination of the dilute regime where 
is the crossover, and (B) vs. where the median value appears as the dotted line. 



with the known interaction potential value. This procedure was applied to all values of r/a, 
creating a data set for each potential that is of the form of <fi m ax vs. r/a; an example data 
set is plotted in Figure [2j3. We should note that in Figure [2j3 it is clear that 4> ma x decreases 
as we increase r/a and it is not simply a single value for a given pair interaction. The 
dependence of <fi m ax on r/a for a fixed e and na can be understood by thinking in terms of 
the "effective" diameter of the particle. That is, if we were to increase the diameter of our 
particle, we would need to increase our system volume to maintain the same (f>. The quantity 
r/a can be considered to be proportional to the "effective" diameter of the particles, thus, 
as we increase r/a we need to also increase our system volume to maintain the same effective 
4>, as we see in Figure [2JB (i.e. increasing system volume is analogous to a decrease <p since 
is calculated using the excluded volume diameter, a, not an effective diameter). As such, 
we calculated the median value of <f) m ax, as this should be representative of the data and not 
biased by the large values of max at small r/a. 

The median values are summarized in Figure [3j The data sets were grouped by e and fitted 
with a power law to guide the eye. We see that as na decreases (i.e. our potential becomes 
longer ranged with respect to particle diameter), 4> m ax,median decreases. Additionally, as e 
increases, we find that 4> ma x,median also decreases, however not as rapidly as the trends with 
na. We also see that for the hard-sphere system, 4> ma x,median = 0.2, and this can be taken 
as an upper limit for all repulsive systems; 4> m ax,median for hard-spheres was calculated using 
the same procedure as the Yukawa systems and corresponds to a state where n = oo, thus 



13 





. 


' 'i ' i 




- 


Hard Sphere Limit, tb .. = 0.2 - 

r ' T max median 










, 


0.1 


_ 


— 


Hard Sphere Limit 


@i< 




: 


D 


£=100 


O^/ : 


c 
CO 
T3 


- 


A 
o 


e=50 

£=3 


-- V 
.-•<£' 


E 




o- „ * • 


X 

CO 




• " *& •' 


c= 




." y 


"^ 0.01 


- 


- x rv 




- 


. -" S ' 






- X 






•* 

X 






X , , , , »• 



0.2 



KO 



FIG. 3. The median value oi(\) max vs. kct showing the approximate crossover between concentrated 
and dilute regimes. Data for each e value are fitted with a power law best fit to guide the eye. 



it is plotted as a bounding line. It is clear that as the range of the potential increases, we 
must sample at increasingly lower volume fractions. 

Figure [3] can be used as a rough guideline to determine the appropriate regime to collect 
data for the calculation of U(r). While this plot depends on a priori knowledge of ko~ and e - 
the values we ultimately wish to calculate - we only need a rough estimate of the screening 
length to use this plot. As previously stated, the median value of (p ma x, median depends 
predominantly on kg which can be roughly estimated with minimal effort by measuring 
at a low value g(r, <p) (e.g. = 0.01 is satisfactory for the systems presented in Figure 
3]). In practice, the values of ma x should be assessed for each system studied using the 
methodology previously described; Figure [3] should be used as a rough guide to approximate 
where this regime occurs, to avoid collecting unnecessary data. 



B. Accuracy of the Extrapolated Potential 



We applied the extrapolation method to our Monte Carlo simulations and found excellent 
agreement between the extrapolated potentials and the pair potentials programmed into the 
simulation; note in these cases we treated the system as monodisperse with respect to particle 
diameter. In Figure [4] we plotted the known potential and the extrapolated potential for an 
example system with e=3 and kg =1.5. We extrapolated W(r,<p) for 11 values of <p ranging 
from 0.002 to 0.03, using the median value of 4>max,median = 0.03 as the maximum cutoff for 
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the regression; we did not take into account the r/a dependence of 4> ma x for this test. We 
see that the extrapolated potential is virtually identical to the known potential, as shown 
in Figure |4j To parameterize the accuracy of the extrapolated potential, we calculated 
the magnitude of the difference between the test and known potentials, normalized by the 
magnitude of the two potentials (represented as vectors), 



Tir / j j \ u test\' i) u known\'i)\ (1Q\ 

M ~ Tr^test itt /._ m , T^max < T t 7~~T7 V id J 



7 j { \Utest{ r i) ~ Uknown\ r i)\ 

/ j \ \Utest\ r i)\ + / jj \Uknown\ r i)\ 

where max is the total number of discrete points considered; in this measure an ideal match 
has M = and the maximum difference has M= 1. We used M rather than the relative 
error since M is always normalized between and 1 for all potentials and number of data- 
points considered. Additionally, by its construction, M is weighted with respect to the total 
magnitude of both vectors we compare, thus we avoid any erroneously large errors associated 
with the difference between relatively small values (e.g. U(r) = 0.05 and U(r) = 0.01 are 
nearly indiscernible on the energy scale we consider here, however their relative error is 4, 
whereas the difference between U(r) = 10 and U(r) = 6 would be clearly evident, yet the 
relative error is much less with a value of 0.67). The difference between the extrapolated 
potential and known pair potential in Figure [4] was calculated to be M=0.023. 

The advantage of using the regression method over a single "dilute" concentration is 
also highlighted in Figure El We plotted W(r, 0) for two different volume fractions, <p 
=0.05 and 0.005. For a value of <fi = 0.05, a relatively low-density system, we capture the 
gross behavior, however the actual values deviate substantially from the known potential. 
W(r, 0.05) significantly undershoots the known potential and any values of e and k extracted 
from this plot would be misleading; for W(r, 0.05), M =0.497. We also plotted W(r, <j)) for 
= 0.005, a very low density sample, and find this to be a reasonable approximation to the 
known potential with M = 0.063. 

Without prior knowledge of the known potential - besides the notion it should be repulsive 
- it would difficult to assess if W(r, 0.05) or W(r, 0.005) are good approximations of the 
pair potential. In order to assess whether the potential of mean force is reasonable, we 
would need to follow a procedure similar to that outlined in reference^ we would need to 
collect additional data at lower values and see if g(r, 0), or alternatively W(r, 0), changed 
substantially as was decreased^. In effect, we would need to collect, analyze, and calculate 
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FIG. 4. For a Yukawa system with e = 3 and kg = 1.5, the pair potential and extrapolated 
potentials are plotted showing excellent agreement. The extrapolated potential was calculated 
from 11 values of <j> ranging from 0.002-0.03. W(r, 4>) is also plotted for 4> = 0.05 and 0.005. M = 
0.023 for the extrapolated potential, M = 0.497 for W(r,0.05), and M = 0.063 for W(r, 0.005) 

additional data that would not necessarily be directly used to calculate the approximation 
of the pair potential. Therein lies a substantial benefit of the extrapolation method. In 
addition to calculating a more accurate approximation of the pair potential, we use nearly 
all data collected to calculate it, thereby increasing our overall confidence in the result; 
this procedure also provides a built-in "check" regarding the appropriateness of the derived 
potential in terms of its dependence on <fr (i.e. if W (r, <fi) does not scale linearly with <p, we 
are not dilute enough). 



In Figure [5j we plotted the extrapolated potentials calculated for e = 3, and ko = (5.0, 
3.0, 0.5, 0.25) to explore the behavior of this method as a function of screening length. As 
in Figure |1J we ignored the r/a dependence of 4> max and used only the median value. The 
extrapolation method worked well in all cases. The average magnitude of the difference 
between the known pair potential (i.e. the potential coded into the simulation) and extrap- 
olated potential for e = 3 and kg = (5.0, 3.0, 0.5, 0.25) is M = (0.028, 0.021, 0.046, 0.087), 
respectively. As the range of the potential increases (i.e ncr becomes smaller), the method 
slightly underestimated the known value of the pair potential at larger r jo values; this was 



most evident for ko = 0.25. As the range of the potential increases, the median of 



or 



any single value, is less representative of the data on a whole, and a more careful treatment 
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of the regression should be applied (i.e. 4> max should be determined for each value of r/a) 
However, the deviations are still well within the expected experimental error. 




FIG. 5. Comparison between the known and extrapolated pair potentials. In all cases, the known 
potential is shown as a solid black line, and the extrapolated potential is shown as circles. In all 
cases, the extrapolated potential matches the known potential showing only minor deviations as 
the range of the potential is increased (i.e. kg decreases). 

In Figure [6] we explored the impact of e on the extrapolation method, plotting the results 
for system where e = 50 and ko = (3.0, 1.0) and e = 100 and ko = (3.0, 1.0). In all cases, 
the extrapolation method matched the known potential well. The average magnitude of 
the difference between the known and extrapolated potentials for e = 50 and ko = (3.0, 
1.0) is M = (0.038, 0.087), respectively and for e = 100 and ko = (3.0, 1.0) M = (0.022, 
0.024), respectively. However, we see that the extrapolation method does not resolve values 
of U(r)/kBT greater than ~10. This is not a failing of the method, but rather a consequence 
of the nature of the strong interactions in these cases. That is, as the strength of the potential 
is increased, the likelihood of two particles coming into very close contact is lowered. 



C. Effect of Polydispersity 



In sections HI A| and IHB we treated our model colloidal particles as ideal and thus 
ignored polydispersity in particle diameter. In practice, experimental systems will contain 
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FIG. 6. Comparison between the known and extrapolated pair potentials. In all cases, the known 
potential is shown as a solid black line, and the extrapolated potential is shown as circles. The 
extrapolated potential matches the known potential well in all cases, however, we can only resolve 
the potential on the order of ~10 fe^T. 

particles with a range of diameters, often distributed with Gaussian behavior, as we observed 
in Figure [I] Polydispersity is denned as, P = (100 %) (standard deviation of particle diam- 
eter) /(average particle diameter). The experimental systems we explored in this work (see 
section II C ) have polydispersity levels of P ~ 4% and this level is not atypical of colloidal 



systems for self-assembly. Polydispersity has been shown to have a large effect on phase be- 
havior, for instance creating a phase separated crystal under high density when P > 8%®l 
However, since we are in the dilute regime, the effects of polydispersity may be muted. In 
our simulations, we initialized our particle diameters based on a Gaussian distribution that 
satisfied the given level of polydispersity we wish to study; our methodology was previously 
described in section III Bt 

We simulated three Yukawa systems with e = 3 and kg = 5, e = 3 and kg = 10, and e 
= 43 and kg = 20, each for P = (0 %, 5 %, 10 %); the extrapolations are shown in Figure 
[7| For e = 3 and kg = 5, polydispersity has a minimal effect on the extrapolation if we 
consider the extrapolated U(r) for r/cr > 1, i.e we do not consider t/g values that are less 
than the separation between two average diameter particles. Specifically, for P = (0,5,10)% 
we find M = (0.028, 0.033, 0.049) respectively when calculating M for r /g > 1. 

Increasing kg to 10 makes the potential shorter ranged and steeper. In this case, poly- 
dispersity appears to have more of an impact; specifically, P =10% undershoots the known 
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potential. However, as we saw for kct— 5, truncating U(r) at r/a > 1 results in nearly 
identical extrapolated pair potentials at all levels of polydispersity. For P=(0,5,10)% we 
find M = (0.030, 0.041, 0.080) respectively when calculating M for r/a > 1. 

For a system with e = 43 and kct = 20, the range of the potential, with respect to the 
surface of the particle, is roughly equal to the case where e = 3 and ko = 10, however it 
increases more rapidly at small r/a. This rapid increase in interaction strength results in a 
system where the diameter of the particle effectively increases (since it is unlikely for particles 
to interact when U(r)/kbT > 10). Hence, the effective interaction is actually shorter-ranged 
than e =3 and kct = 10 (i.e. the interaction can be thought of as shorter-ranged and radially 
shifted outwards from the surface of the particle). This is particularly evident as even the 
P =0% system has trouble resolving the known pair potential well at small r/a. The effect 
of polydispersity seems more pronounced than the other systems, undershooting the known 
potential in both the 5% and 10% cases. In this case, simply truncating our extrapolation 
at r/a > 1 is not sufficient to minimize the impact of polydispersity. For P=(0,5,10)% we 
have M = (0.302, 0.47, 0.59) respectively when calculating M for r/a > 1. However, we 
note that for r/a > 1.125, we are able to resolve the potential reasonably well for all levels 
of polydispersity. 



From this we can conclude that as the range of the potential decreases, the effect of 
polydispersity increases. This can be understood as arising from the fact that as the potential 
acts over shorter distances, the likelihood of close contact increases. In a polydisperse system, 
where a subset of particles have diameters smaller than the average, we will observe a greater 
number of particles at small separations of g(r, <fi) within the excluded volume region. Thus, 
we can expect that relatively short-ranged potentials (or those that are effectively short- 
ranged) will show a stronger effect of polydispersity and our measurements of e will be less 
accurate than k. As a general rule, we should exclude any values of the derived potential that 
occur for r/a < 1 and apply extra scrutiny to any values at small r/a, especially in systems 
where the derived potential appears short-ranged. Additionally, any potential estimation 
method that relies on g(r,<f)) will suffer from this same issues regarding polydispersity and 
very short-ranged potentials. 
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FIG. 7. Comparison of extrapolated potentials with different levels of polydispersity. The grey 
line denotes the excluded volume region, which corresponds the minimum separation between two 
average sized particles. 

IV. EXPERIMENTAL RESULTS 



A. Analysis of Pure DOP System 



To test the validity of the potential derivation method in experiment, we first considered 
PMMA/PHSA particles in pure DOR We collected data for a range of = (0.005, 0.007, 



0.008, 0.012, 0.018, 0.022, 0.042), as summarized in table IV A For the three lowest values 
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of — (0.005, 0.007, 0.008), a large number of samples volumes were used for statistical 
purposes (see appendix regarding error scaling in g(r)). Figure [8] plots the average g(r, 0) for 
= (0.005 - 0.022) where the bin shell size Ar = 0.1/xm ( 10% of the particle diameter) to 
balance the competing constraints of signal to noise and spatial resolution of the potential. 
At these charge and solvent conditions, we observed no particles in the inner most shells 
of g(r,4>). This result is consistent with repulsive pair interactions, as reported for similar 
systems^^. Error bars represent the standard error of the mean of g(r, 0). 
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FIG. 8. g(r, <f>) for PMMA/PHSA particle system in pure DOP for = (0.005, 0.007, 0.008, 0.012, 
0.018, 0.022) respectively from the bottom. 



Figures[9^-C plot —ln[g(r, 0)] (i.e. W(r, 0)) as a function of for 0.005 < < 0.042 for r 
= (1.75, 1.85, 1.95) /zm for the pure DOP system. As we saw in the simulations (see Figure 
2l), there is a clear region that displays linear behavior; the crossover between "dilute" and 
"concentrated" appears at ~ 0.0125 which is nearly identical to the value of 4> ma x,median ~ 
0.015 approximated from Figure [3] for e =27 and kg =1.1 from table |ITCT 

The pair potential calculated using the extrapolation method for PMMA/PHSA in the 



pure DOP solvent is plotted Figure 10 U(r)/kBT is qualitatively consistent with long-range 
repulsive interactions, as reported for similar systems^. Error bars plotted are standard 
error showing the data is precise. The theoretical description of the potential (derived 
from electrokinetics and the one parameter fit of the unknown solvent conductivity) is also 



plotted in Figure 10 as reported in table II C 1 These two potentials are in agreement; the 
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N total 


Total # of 
image volumes 


0.005 


7216 


36 


0.007 


10911 


26 


0.008 


5772 


12 


0.012 


4686 


6 


0.018 


6842 


6 


0.022 


8279 


6 


0.042 


15997 


6 



TABLE II. Pure DOP solution. N total refers to the total number of particles from the multiple 
image volumes used in the calculation of g(r, (j>). 

average magnitude of the deviation between the extrapolated potential and the theoretical 
prediction is M = 0.2. Note that we expect the value of M to be greater for the experimental 
system than simulation due to experimental noise; e.g. M— 0.15 when comparing the raw 
data to a best fit of the form of the Yukawa potential. We additionally plot W(r, 0.022) 
and W(r, 0.008). As we saw in our simulations, W(r,<p) calculated from even relatively 



low <p values may drastically underestimate the potential (see Figure 10), and the only way 
to determine this would be to collect additional datasets. In Figure 10 we plot W(r, 0.022) 
where M = 0.65, noting that even this relatively low volume fraction provides a poor estimate 
of the potential. We also plot W(r, 0.008) which matches both the extrapolated and theory 
potential well, with M = 0.18. 

B. Effect of TBAC Addition 

To assess the extent to which the method can resolve differences in pair potential inter- 
actions that are relevant for self-assembly, we performed experiments in which the range 
of repulsive interactions are reduced by the addition of electrolyte. We followed the same 
procedure previously outlined. We collected data at <p = (0.012, 0.019, 0.024, 0.037, 0.049) 
for the 10/iM TBAC system and at = (0.012, 0.018, 0.031, 0.039, 0.052) for the 2mM 
TBAC system. For both TBAC concentrations, we observed that the onset of finite g(r,<f>) 
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FIG. 9. W(r,(f>) vs. <fi for radial positions of (A) r = 1.75 ^m, (B) r = 1.85 /im, and (C) r = 1.95 
/im for the pure DOP system. The linear, "dilute" regime appears for eft < 0.0125. 



is shifted to smaller radial distances (see appendix for more detail). This shift qualitatively 
indicates that TBAC addition moderates repulsive interactions in the system. This shift also 
affects the limit of linearity. For example, for the two cases of added salt ( [TBAC] = 10 
/jjM and 2 mM) we found that linearity was maintained out to <fi ~ 0.05 which corresponded 
to the highest volume fraction studied. This increase in the limit of linearity is consistent 
with the simulations, for which <p max, median ~ 0.1 from Figure [3j That is, compared to the 
case of pure DOP, where linearity is found only for < 0.0125, these results coincide with 
an expected shift due to an increase in the screening of the repulsive interactions shown in 
simulation. 



The extrapolated pair potentials for the systems with TBAC are plotted in Figure 11 
Visually, we find good agreement between the theoretical potential calculated from elec- 
trokinetics for 10/iM system where M = 0.35, as shown in Figure 11 \ . For [TBAC] = 2 



mM, we find that the extrapolated potential does not provide as close a match as the pure 
DOP and 10/iM systems when compared to electrokinetic-based potential; for [TBAC] = 2 
mM, M = 0.84. The 2mM TBAC system is nearly consistent with simple excluded volume 
interaction; the transition to the repulsive portion of the potential is abrupt and very close 
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FIG. 10. U(r) for PMMA/PHSA particle system in pure DOP. Screened Coulomb potential cal- 
culated from electrokinetics/fitting is shown as a solid line, the extrapolated potential is shown as 
open circles. W(r, 0.022) is plotted as filled circles, and W(r, 0.008) is plotted with x symbols. The 
gray line denotes the excluded volume region appearing at r < 1 ^m. 



to the measured diameter of the colloid, with a Debye length of 49 nm. Thus, screened 
electrostatic interactions extend no further than about 5% of the particle diameter, not too 
different from the spread in the particle size distribution due to polydispersity which was 
calculated to be 4%. The effect of polydispersity on this 2mM system is consistent with what 
was observed in simulation, as discussed previously in Figure [7j as the range of the potential 
decreases, polydispersity plays a stronger role and the observed e value decreases. However, 
if we truncate the extrapolated potential at r ~ 1.25/im, ignoring any interactions within 
the excluded volume region and close to the surface of an average particle, we arrive at a 
more satisfactory description of U(r). We should also note that, as seen in the simulations, 
we do not resolve U(r)/kBT > 10 for any of the systems studied; however, even given this 
constraint we can still easily resolve the shape of the potential. 
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FIG. 11. U(r) for PMMA/PHSA particle system for (A) 10/iM TBAC, and (B) 2mM TBAC. 
Screened Coulomb potentials calculated from electrokinetics/fitting are shown as solid lines, ex- 
trapolated potentials are shown as symbols. The gray line denotes the excluded volume region 
appearing at r < 1 p. 

V. CONCLUSIONS 



We have presented a method to calculate the pair potential, U(r), between colloidal 
particles from microscopy. This method relies on extrapolating the potential of mean force, 
W(r,(j)), to = 0. We have shown using MC simulation that this method produces near 
perfect results for ideal monodisperse colloids. We have used simulation to explore the 
impact of surface charge and screening length on this method, providing general guidelines 
for the use of this method over a wide range of parameters. We have also demonstrated 
that low levels of polydispersity have only a small impact on the accuracy of the calculation 
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of U(r) for longer-ranged potentials, however the effect is more pronounced as the range 
of the potential decreases. Further, we have applied this method to experimental colloidal 
particles using the guidelines established from simulation. We found close agreement of 
the results of the extrapolation method with theoretical screened Coulombic potentials. 
As was noted in simulation, the effect of polydispersity becomes stronger for short ranged 
potentials where the range of the potential is approximately equal to the spread in particle 
diameter. Moreover, we see that this methodology is well suited to determine the potential 
for a range of systems, including those that are refractive index matched, particularly for 
systems where the screening length is on the order of the particle diameter and many body 
effects are minimal. We anticipate that this methodology can also be used to quantify 
U(r)/kbT for attractive systems, assuming the attraction is not so strong as to induce 
gelation. Additionally, this method can be used in concert with electrokinetic measurements 
as a means to double-check their validity or used in their place when it is not possible to 
measure key parameters (e.g. conductivity). 
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VII. APPENDIX 

A. Addition of TBAC to the DOP solvent 



In section IV B we explored the addition of TBAC to the DOP solvent. The total num- 



ber of image volumes/particles collected are summarized in table VII A for the 10 /zM solu 



tion and|VIlA for the 2mM solution. The g(r,(j)) data used to calculate the extrapolated 



U(r)/kbT (extrapolations were previously shown in Figure 10) are plotted for in Figure 12 



<t> N tota i Total # of 
image volumes 

0.012 6646 10 

0.019 8331 8 

0.024 10434 8 

0.037 12183 8 

0.049 16116 6 

TABLE III. [TBAC] = 10 /iM. 



B. Error scaling in the radial distribution function 

Essential to accurately determining the pair potential is generating g(r, <fi) with sufficient 
accuracy. g{r,<f)) is a measure of the probability of finding a particle at a distance r away 
from a given reference particle. This distribution function is determined by calculating 
the distance between all particles pairs and binning them into a histogram; the histogram 
is then normalized with respect to an ideal gas, where particle histograms are completely 
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<P 


N total 


Total # of 
image volumes 


0.012 


10568 


14 


0.019 


13508 


12 


0.031 


10196 


6 


0.037 


12453 


6 


0.049 


15129 


6 



TABLE IV. [TBAC] = 2 mM. 




r (urn) 



r(nm) 



FIG. 12. g(r,(f>) plots for [TBAC]= 10 /jM (left) and [TBAC] = 2 mM (right). Data corresponds 
to <f> = (0.012,0.019,0.024, 0.037, 0.049) [TBAC] = 10 /M (left)and <j) = (0.012,0.019,0.031, 0.037, 
0.049) for [TBAC] = 2 mM (right) respectively from the bottom. 



uncorrelated. Due to its construction, g(r, <fi) is a spherically averaged measure; therefore, we 
average over more particle pairs as we consider larger particle separations. In other words, 
our statistics increase as r increases. This may be problematic since the pair potential we 
wish to determine acts the strongest at smaller separations where we have the least statistics. 
Additionally, since we must compute g(r, <fi) at low values of 0, it is unlikely to have many 
particles interacting at small values of r if particles are well distributed in the system (i.e. 
particles do not aggregate); unlike in a dense system, a single configuration snapshot will 
not be sufficient to accurately generate g(r,(f)). As such, to determine U(r) accurately, we 
must establish when we have appropriately gathered sufficient statistics to be confident in 
9(r,(f>). 

We considered the error in g(r, <fi) at a given separation of r to be proportional to l/V-W&iru 
where Ni> in is the number of particles in the bin. By this construction, our maximum error 
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FIG. 13. Error scaling in g(r,cj)) for a combination of Yukawa systems with e = 3 and 
k = (3.0,1.5,0.5,0.25), over a range of volume fractions, 4> = (0.001,0.005,0.01,0.02, and 
Ar = (0.05,0.15,0.25). 



is 1/a/2 = 0.707 and our minimum error will approach 0. There are three main factors 
that impact this value: the error decreases as we increase the total number of particles 
considered, N tota i (i.e. the sum of all particles in all the samples used to construct g(r, </>)); 
the error decreases as we increase our volume fraction, <b\ and the error decreases as we 



increase the width of the bin, Ar. In Figure 13 we plot 1/y/Nun vs. (N to tal) (</>){ Ar) for 
a combination of Yukawa systems; here we have chosen to plot the data for a bin close to 
the surface of the particle, specifically r jo = 1.125, as this should be characteristic of the 
maximum error in the system. We simulated four different potentials; we fixed e = 3 and 
varied the inverse screening length, specifically k = (3.0, 1.5, 0.5, 0.25). For each potential 
we performed simulations over a range of volume fractions, = (0.001, 0.005, 0.01, 0.02) 
and then analyzed each of these simulations using three bin sizes, Ar = (0.05, 0.15, 0.25). 



We plotted the data from these four potentials in Figure 13 finding that the error roughly 



scales as a power law with exponent -0.57. The inset of Figure 13 plots the data on a 
standard axis; it is clear from this plot that our error is most rapidly decreasing in the range 
of 1 < (N tota i)(<p)(Ar) < 10. At values of (iV toia ;)(0)(Ar) > 10, our error decreases very 
slowly. Thus for the best statistics in the g(r,4>), a rough guideline should be to perform 
enough samples such that (iV to4a j) (</>)( Ar) > 10. For our simulations Ar = 0.1 and (N to tai = 
1000000, thus, even for our lowest volume fraction of d>= 0.001, this condition is met. 
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Similarly, Figure 14 summarizes the relationship between the precision of the measured 
g(r, <fi) and properties of the experiments, including the number of particles, the volume 
fraction of the specimen, and the bin size of each point in g(r,(j>). We see that the error 
in g(r,(j)) is a weak function of the combination (iV tclta ;) (</>)( Ar), scaling as power law with 



exponent -0.79. Figure [14] can be applied to determine experimental conditions for pair po- 
tential characterization. For example, for a relative error of ~10% in the radial distribution, 
Figure 14 requires (N tota i)(<j))(Ar) ~ 7. For a volume fraction 0.01 with 1 /zm particles in 



which Ar = 0.1 /xm, this relative error would require on the order of 7000 particles. 




(N)fo)(Ar) 



FIG. 14. Error scaling in g(r, 4>) as a function of number of particles, volume fraction, and bin size 
for pure DOP, no salt case. Data taken from four different radial positions, r = (2.05, 2.15, 2.35, 
2.55, 2.75, 3.05) /un. Line drawn is power law fit with exponent -0.79. 
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